#8305211218-吴滔
#代码改自课件
library(dplyr)
library(ggplot2)
library(Seurat)
##单细胞测序####
#质量控制和数据过滤
setwd("D:\\R_DATA")
pbmc.data <- Read10X(data.dir =  "D:\\R_DATA\\filtered_gene_bc_matrices\\hg19")#读取文件，默认读取基因为第二列
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
#创建Seurat类对象，基因必须要在3个及以上细胞中表达，细胞必须要有至少200个基因表达）
pbmc[['percent.mt']] <- PercentageFeatureSet(pbmc, pattern = "^MT-|^mt-")#创建percent.mt列用于表征线粒体基因表达
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#用小提琴图看下表达情况，确定过滤阈值
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
#数据处理，过滤条件:单个细胞中至少检测到200个基因且不多于2500个基因，线粒体基因比例不超过5%
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#数据标准化处理，标度因子为10000，转换为log1p（去除测序深度带来的影响）
all.genes <- rownames(pbmc)#提取基因列表
pbmc <- ScaleData(pbmc, features = all.genes)#数据缩放
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)#寻找高差异表达基因，方法为方差稳定变化（默认值）
pbmc <- ScaleData(pbmc)#数据中心化

#额外操作
#储存前10位高变基因
top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)#高变基因散点图
plot2 <- LabelPoints(plot=plot1, 
                     points = top10,
                     repel=TRUE)#top10加上基因名标签
CombinePlots(plots = list(plot1, plot2), ncol =2)#结合到一张图中

#主成分分析
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
ElbowPlot(pbmc)#查看方差骤降图
DimPlot(pbmc, reduction = "pca")#绘制降维散点图
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")#可视化降维贡献


#细胞聚类和绘制UMAP图####
pbmc <- FindNeighbors(pbmc, dims = 1:10)#选择前10维度
pbmc <- FindClusters(pbmc, resolution = 0.5)
pbmc <- RunUMAP(pbmc, dims = 1:10)
my_umap<-UMAPPlot(pbmc)

#分析差异表达
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
#正向标记基因-25%的细胞检出率-fc阈值
top2 <- pbmc.markers %>% group_by(cluster) %>% slice_max(n = 2, order_by = avg_log2FC)
#发掘每个cluster top2基因

#气泡图-feature plot和VlnPolt图####
p10<-DotPlot(object =pbmc ,assay = 'RNA',features = top2$gene)+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
#绘制气泡图，X轴标签旋转45度，右缘对齐
mytop<-top2[c(1,3,5,6),]#挑选几个marker基因
myfp<-FeaturePlot(pbmc, features = mytop$gene)#FeaturePlot图
myvp<-VlnPlot(pbmc, features = mytop$gene)#VlnPlot图


#统一输出pdf文件
pdf("UMAP.pdf",width=8,height = 6)
print(my_umap)
dev.off()
pdf("Cluster_DotPlot.pdf",width=8,height = 6)
print(p10)
dev.off()
pdf("FeaturePlot.pdf",width=8,height = 6)
print(myfp)
dev.off()
pdf("VlnPlot.pdf",width=8,height = 6)
print(myvp)
dev.off()
